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We present the self-consistent implementation of current-dependent (hybrid) meta generalized gradient 
approximation (mGGA) density functionals using London atomic orbitals. A previously proposed generalized 
kinetic energy density is utilized to implement mGGAs in the framework of Kohn-Sham current density- 
functional theory (KS-CDFT). A unique feature of the non-perturbative implementation of these functionals is 
the ability to seamlessly explore a wide range of magnetic fields up to 1 a.u. (~ 235000T) in strength. CDFT 
functionals based on the TPSS and B98 forms are investigated and their performance is assessed by comparison 
with accurate CCSD(T) data. In the weak field regime magnetic properties such as magnetizabilities and 
NMR shielding constants show modest but systematic improvements over GGA functionals. However, in 
strong field regime the mGGA based forms lead to a significantly improved description of the recently 
proposed perpendicular paramagnetic bonding mechanism, comparing well with CCSD(T) data. In contrast 
to functionals based on the vorticity these forms are found to be numerically stable and their accuracy at high 
field suggests the extension of mGGAs to CDFT via the generalized kinetic energy density should provide a 
useful starting point for further development of CDFT approximations. 
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I. INTRODUCTION 

The foundations of current density-functional theory 
(CDFT) and its Kohn-Sham (KS) implementation were 
established in the late 1980’s with the seminal works of 
Vignale, Rasolt and Geldart 1-3 , where it was recognised 
that the exchange-correlation functionals must depend 
not only on the electronic density, p , but also the param¬ 
agnetic current density j p in the presence of an electro¬ 
magnetic field. Since these early works a large number of 
theoretical investigations of CDFT have been presented. 
The foundations of CDFT have sometimes been viewed 
as controversial. Most recently, Pan and Sahni 4-6 sug¬ 
gested that the physical current density j, rather than 
j p , aught to be the fundamental variable in a CDFT and 
attempted to establish a Hohenberg-Kohn like theorem 
for this physically appealing alternative choice of variable. 
Unfortunately, the derivation of this theorem has been 
shown to be in error' ,8 . Furthermore, the work of Tellgren 
et al.‘ showed how CDFT may be brought into Lieb’s 
convex-conjugate formulation of DFT 9 , further strength¬ 
ening its foundations and lending key insight into the 
more complex relationship between the key densities and 
potentials in the theory. In particular, it is highlighted 
that the lack of a Hohenberg-Kohn theorem is not an 
impediment to a viable CDFT. Recent theoretical works 
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by Lieb and Schrader 10 and Tellgren et al . n have also 
addressed the issue of ./V-representability in CDFT. 

Despite the theoretical progress in CDFT very few prac¬ 
tical implementations of theory have been presented. Most 
practical studies have either presented calculations based 
on fixed densities (typically computed at the Hartree-Fock 
or standard Kohn-Sham level), or have attempted to in¬ 
clude CDFT contributions in linear-response calculations. 
In the context of response theory, implementations have 
been presented for magnetic properties by Lee et al. 12 
and for excitation energies at the meta generalized gradi¬ 
ent approximation (mGGA) level by Bates and Furche 13 . 
Very few fully self-consistent implementations of CDFT 
capable of treating systems beyond the linear-response 
regime have been presented, in fact we are only aware 
of the work by Pittalis et a/. 14,15 in the context of the 
optimized effective potential method, Zhu and Trickey 16 
for atomic systems and our own implementation 1 ' for 
general atomic and molecular species. 

A number of challenges arise when implementing quan¬ 
tum chemical methods for molecules in magnetic fields, 
the London program 18 has been specifically designed to 
address these and is utilized throughout the present work. 
In particular, London atomic orbitals 19-21 are employed 
to ensure gauge origin invariant results. For CDFT addi¬ 
tional challenges arise since new forms are required for the 
exchange-correlation functional. Relatively few practical 
forms for CDFT functionals have been suggested in the 
literature. 

In the present work we will examine the use of mGGAs 
and hybrid mGGAs for the exchange-correlation energy 
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in the presence of magnetic fields. Functionals in this 
class depend on the orbital dependent kinetic energy den¬ 
sity in the absence of a magnetic field. However, as has 
been noted in the literature 13,22,23 , this key quantity is 
not gauge invariant and so some modification is required 
for use in a magnetic field. One approach is to replace 
the kinetic energy density by a generalized form including 
the paramagnetic current density. This quantity natu¬ 
rally arises in the expansion of the spherically averaged 
exchange hole, as derived by Dobson 24 . Becke has already 
suggested the use of this approach to produce a current de¬ 
pendent generalisation of the Becke-Roussel 25 functional. 
Recently, Bates and Furche 13 have also explored a simi¬ 
lar generalization of the Tao-Perdew-Staroverov-Scuseria 
(TPSS) functional 26 to calculate excitation energies via 
response theory. 

We will consider current dependent extensions of the 
B98 27 , TPSS 26 and TPSS(h) 28 functionals. The use of 
a modified current-dependent kineitc energy density is 
denoted by a prefix c throughout the remainder of this 
work. The non-perturbative nature of the implementa¬ 
tion in the London program will allow for testing of 
these functionals in both weak and strong field regimes. 
The availability of accurate ab initio methodologies in 
the London program provides a unique opportunity for 
the assessment and testing of CDFT functionals at field 
strengths upto 1 a.u. In the present work we make use 
of the recent implementation by Stopkowicz et al. 29 of 
coupled-cluster (CC) methods with single, double and per¬ 
turbative triple excitations [CCSD(T)] for benchmarking 
the CDFT approximations. 

In Section II we review the simple generalization of 
mGGA functionals to the CDFT framework due to Dob¬ 
son and Becke, details specific to the functionals con¬ 
sidered in this work are collected in the appendix. In 
Section III we outline some computational details of our 
calculations. Section IV summarizes our findings, assess¬ 
ing the quality of the CDFT approximations by com¬ 
parison with CCSD(T) data; in Section IV A we explore 
the performance of mGGA functionals for calculating 
molecular properties in the weak held regime accessible 
via linear response theory; in Section IV B the high held 
regime is explored by considering the recently proposed 
perpendicular paramagnetic bonding. The interpretation 
of this bonding mechanism in the KS-CDFT framework 
is discussed in Section V. Finally, concluding remarks and 
directions for future work are given in Section VI. 


II. THEORY 

In this work we consider the calculation of energies and 
molecular properties in the presence of a static uniform 
external magnetic held, B, which may be represented in 
terms of the vector potential 


where Rg is an arbitrary gauge origin. The London 
program makes use of London atomic orbitals 19-21 to 
ensure that computed energies and molecular properties 
are invariant with respect to choice of the gauge origin. 
These basis functions take the form 


w M (rK-,B,R G ) = exp 


-B x (Rq — Rk) • r 


X^k) 


( 2 ) 

where is a standard Gaussian-type orbital centred 

at position Ra-. These perturbation-dependent basis 
functions are used to expand the Kohn-Sham molecular 
orbitals. 

The Kohn-Sham approach to density-functional the¬ 
ory can be extended to CDFT by searching for a non¬ 
interacting system of electrons with the same charge 
and current densities as the physical interacting system: 
(/9 s (r),j PjS (r)) = (p(r),jp(r)). The KS equations can be 


written as 



-{p,A s } + u s 


[V x A s ] 


(Pp — SpPp (3) 


where the KS potentials (u B , A s ) are defined as 

Us = U e xt + Uj + U X c + —A s (4) 

and 


A s — A ext 


+ A X c 


( 5 ) 


The first three terms in the scalar potential u s represent 
the external (ext), Coulomb (J) and exchange-correlation 
(xc) potentials defined as the functional derivative of the 
respective energies with respect to the density, as in the 
usual KS approach. The final contribution to the scalar 
potential arises from the non-interacting vector potential 
A s . In addition to the vector potential due to the external 
field A ex t as defined in Eq. (1) the KS vector potential 
contains an exchange-correlation contribution defined as 


_ ^-Exc [p, jp] 

<5j 


A central question that immediately arises in CDFT is 
how the exchange-correlation functional must be modi¬ 
fied to include current effects. Whilst the paramagnetic 
current density is a valid quantity on which to base the 
universal density functional it can also be shown that the 
exchange-correlation component must be independently 
gauge-invariant 2 . This places a significant constraint on 
the manner in which this quantity may enter any approx¬ 
imate CDFT functional. In contrast to standard DFT, 
relatively few CDFT functionals have been proposed. The 
majority of these are based on the the vorticity 


with 


i/(r) = V x 



( 7 ) 


A(r) = ^B x (r-Rd), 



SE xc [p,v] 
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as proposed by Vignale, Rasolt and Geldart (VRG) 3 . The 
original VRG form for the exchange-correlation energy 
was parameterized using Monte Carlo simulations of the 
high density limit 1 . A number of re-parameterisations 
for this form have been suggested based on accurate cal¬ 
culations in the high-density regime 12,30-32 . Higuchi and 
Higuchi 33 (HH) have also presented a vorticity dependent 
form, derived to obey known exact relations for the CDFT 
exchange-correlation functional. 

Whilst the vorticity is a theoretically convenient choice 
to ensure the gauge invariance of the exchange-correlation 
energy it has been observed that in practical self-consistent 
calculations it can lead to significant numerical stability 
issues 16,1 ' . How severe these issues are depends on the 
exact parameterization of the the functional form, how¬ 
ever, in all cases some degree of numerical regularization 
is required to ensure that the self-consistent field solution 
of the Kohn-Sham equations can be obtained. Further¬ 
more, molecular properties computed by such calculations 
display an un-acceptably strong dependence on the reg¬ 
ularization parameters - with no obvious convergence 
towards a single value. Clearly this raises questions as 
to how appropriate such forms are for use in quantum 
chemical calculations. 

Most practical mGGAs make use of the kinetic energy 
density 


occ 

Tt = V ip* a ■ V ipi<j , (9) 

i 


in their construction, where (pi are the occupied KS or¬ 
bitals and a is the electron spin index. This term is gauge 
dependent and as a result an unmodified meta-GGA type 
functional form cannot be used to describe a system with 
a non-zero magnetic vector potential. To resolve this issue 
the gauge independence of the exchange-correlation func¬ 
tional must be restored. A natural modification, which 
can be applied to any mGGA dependent on the kinetic 
energy density r(r), arises in the work of Dobson 24,34 who 
generalized the expansion of the exchange-hole to include 
the case of non-zero current densities. 

The spherically averaged exchange hole at zero field can 
be modelled using a Taylor expansion 35 and is commonly 
considered 25,27,36 up to the quadratic term 


Qa 


1 

6 


V 2 /? ct — 2t ct + 


(V/v) 2 ~ 

2 Pa 


( 10 ) 


This expansion can be generalised to non-zero field and 
the curvature term becomes 24 ’ 25,34 


Qa = l 
6 


V 2 /v - 2 T a 


(Vp.) 2 , 2 |j 


per | 


2p<7 PcT 

where j pcr is the paramagnetic current density 


( 11 ) 


■W = - \ E . (i2) 


Comparing Eqs. (10) and (11) it is possible to identify a 
correction to the conventional r(r) that is gauge invariant 
and may be utilized in mGGA functionals, 


Ta = T a 


IJpcxI 2 

Pa 


(13) 


The use of Eq. (13) has been put forward many times 
in the literature. Becke suggested its use in the Becke- 
Roussel model 25 to generate a current dependent analogue 
of this functional. He also suggested that this quantity 
could be used to define a current dependent inhomogene¬ 
ity parameter in the more empirical B98 functional 2 '. It 
has also been suggested for use to generalize the TPSS 
functional 26 by Tao 23 . Recently Bates and Furche 13 con¬ 
sidered the application of the resulting cTPSS functional 
in the calculation of excitation energies via response the¬ 
ory. In the present work we consider the application of 
mGGA functionals with this modification to calculation 
magnetic properties in the weak and strong field regimes 
in a non-perturbative manner. In particular we consider 
three functional forms, cB98, cTPSS and cTPSS(h), where 
the prefix c denotes the use of the modified f a in Eq. (13). 
The Appendix gives some details of these functional forms 
to show how these modifications enter. 


III. COMPUTATIONAL DETAILS 

Unless otherwise indicated all calculations in this work 
use the aug-cc-pCVTZ basis set 37,38 . All DFT calcula¬ 
tions have been performed using the London 18 program. 
This code utilizes the XCFun library 39 for the evaluation 
of the density functionals and their derivatives. The mod¬ 
ifications of Eq. (13) and the functionals cB98 and cTPSS 
have been added to the XCFun library. In addition we 
investigate the use of a hybrid form of cTPSS, denoted 
cTPSS(h), based on the TPSS(h) functional of Ref. 40. 
The quality of the CDFT functionals cB98, cTPSS and 
cTPSS(h) is assessed by comparison with CCSD(T) data. 
For comparison Hartree-Fock (HF), local density approx¬ 
imation (LDA), Perdew-Burke-Ernzerhof (PBE) 41 , and 
Keal-Tozer-3 (KT3) 42 density-functional results are also 
presented. The latter is of particular interest since it is 
specifically designed for the calculation of nuclear mag¬ 
netic resonance (NMR) shielding constants. 

The performance of these approximations will be con¬ 
sidered in two regimes; the weak field regime accessible 
by linear response calculations and the strong field regime 
only accessible via non-perturbative calculations. In the 
weak field regime we will consider the calculation of molec¬ 
ular magentizabilities and NMR shielding constants for 
the 26 small molecules in Table I. Errors for these quan¬ 
tities are presented relative to the benchmark data of 
Ref. 43. Results are also compared with those includ¬ 
ing corrections from the Tao-Perdew parameterization 31 
of the Vignale-Rasolt-Geldart functional 2 , taken from 
Ref. 17. In the strong field regime we consider the pre¬ 
diction of perpendicular paramagnetic bonding 44 in a 
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TABLE I. The test set of molecules for which accurate bench¬ 
mark CCSD(T) data from Ref. 43 was available. 

HF CO No H 2 0 HCN HOF LiH 

NH 3 H 2 CO CH 4 C 2 H 4 A1F CH 3 F C 3 H 4 
FCCH FCN H 2 S HOP HFCO H 2 C 2 0 LiF 
n 2 o OCS h 4 c 2 o PN so 2 


field strength of 1 a.u. perpendicular to the internuclear 
axes of H 2 , He 2 , HeNe and Ne 2 . These non-perturbative 
calculations are assessed against CCSD(T) results com¬ 
puted using the implementation of Ref. 29 in the London 
program. 


IV. RESULTS 

A. The weak field regime: magnetic properties 

We commence by considering the molecular magneti- 
zabilities and NMR shielding constants of the 26 small 
molecules in Table I. The magnetizability tensor elements, 
are defined as 


£«,/3 


d 2 E{ B) 


dB a dBp 


B—0 


(14) 


where a and /3 label cartesian components of the tensor 
and magnetic field. The NMR shielding tensor for a given 
nucleus K is defined by 




d 2 E( B,M k ) 
dB a dM K ,p 


B—O.Mk—0 


(15) 


where Mjf is the nuclear magnetic moment of nucleus 
K. These properties can be accessed non-perturbatively 
in the London program by explicit calculation of the 
energy in the presence of the perturbing fields. Details 
of this procedure are given in Ref. 17, here we compute 
values for each method considered in the same manner, 
facilitating a comparison with previous results. 

Given that for many density-functional approximations 
these singlet second order magnetic response properties 
can be accessed by standard linear response methods in a 
variety of programs, this approach may seem cumbersome. 
However, it should be noted that the implementation of 
the new CDFT approaches in this framework is much more 
straightforward, requiring only an implementation of the 
functional and the derivatives required for construction 
of the KS matrix. More importantly, as we will see in 
Section IV B, this non-perturbative approach allows us to 
seamlessly explore the behaviour of new approximations 
in much stronger fields - inaccessible via linear response 
theory. This means that London provides a powerful 
testbed for new CDFT functionals. 

To quantify the accuracy of the DFT approaches for the 
calculation of these properties we compare our results with 
the CCSD(T) benchmark values of Ref. 43. Specifically, 


TABLE II. The mean error (ME), mean absolute error (MAE), 
and standard deviation (SD) of magnetic properties relative 
to the CCSD(T) benchmark data of Refs. 43,45 . See also Figs. 
1 and 2. 



Magnetizability 

(io -30 jt- 2 ) 
ME MAE SD 

NMR Shielding 
(ppm) 

ME MAE SD 

HF 

-3.06 

6.35 

7.29 

-15.18 

21.40 

40.97 

LDA 

5.01 

9.18 

10.86 

-24.81 

24.85 

30.00 

PBE 

6.75 

8.81 

9.03 

-19.66 

19.78 

21.46 

PBE+VRG (TP) 

7.85 

9.61 

9.44 

-20.33 

20.46 

22.43 

KT3 

8.18 

8.95 

7.83 

-6.53 

8.94 

13.13 

KT3+VRG(TP) 

9.18 

9.83 

8.19 

-7.45 

9.18 

13.37 

B97-2 

5.46 

5.84 

5.96 

-16.34 

16.48 

20.60 

cB98 

0.52 

4.84 

6.58 

-12.44 

12.66 

17.79 

cTPSS 

7.13 

7.51 

6.76 

-14.14 

14.35 

15.61 

cTPSS(h) 

6.41 

6.51 

6.00 

-14.33 

14.52 

16.68 


we use the values at the CCSD(T)/aug-cc-pCV[TQ]Z level 
- which have been extrapolated to the basis set limit using 
the procedure of Refs. 43 and 45. 

The errors in the calculated magnetizabilities and NMR 
shielding constants are summarised in Table II and pre¬ 
sented graphically in Figures 1 and 2 as box-whisker plots. 
In these plots individual points represent the errors for 
each system relative to the reference values, the upper and 
lower fences of the whiskers denote the maximum positive 
and negative errors respectively and the coloured boxes 
enclose errors between the 25% and 75% quantiles. Mean 
and median errors are marked in the plots by horizontal 
black and white lines, respectively. Grey diamonds are 
used to represent the confidence intervals. 

For the molecular magnetizabilities it is clear from the 
error measures in Table II and their representation in 
Figure 1 that none of the functionals offers high accuracy. 
The GGA functionals PBE and KT3 in particular do not 
offer significant improvements over LDA. Whilst their 
minimum and maximum errors are slightly improved, the 
mean errors actually deteriorate. Similar observations 
were made in Ref. 45 for these type of functionals. The 
B97-2 functional gives slightly reduced errors and this is 
consistent with previous conclusions 45 that for magnetiz¬ 
abilities the inclusion of HF exchange may be beneficial. 
At the GGA level the underlying functionals are already 
gauge invariant but do not depend explicitly on the param¬ 
agnetic current density. To introduce this dependence the 
VRG functional may be added. In our earlier work 1 ' we 
found that this correction can be numerically problematic 
and that the most stable parameterization of this func¬ 
tional to date is that put forward by Tao and Perdew 31 , 
denoted VRG(TP). For comparison we include here the 
PBE+VRG(TP) and KT3+VRG(TP) results. It is clear 
that the inclusion of the VRG (TP) correction actually 
worsens the agreement of the results with CCSD(T). At 
the rnGGA level the inclusion of current is mandatory 
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FIG. 1. Box-whisker plot of errors in (C)DFT molecular 
magnetizabilities relative to the extrapolated CCSD(T) values 
of Ref. 45. All calculations use the aug-cc-pCVTZ basis set. 
See the text for further details. 

to ensure the exchange-correlation evaluation is gauge 
independent. 

Here we investigate the cB98, cTPSS and cTPSS(h) 
functionals, the TPSS and TPSS(h) forms are similar in 
performance to PBE - offering marginal improvements on 
some error measures, with TPSS(h) performing slightly 
better than TPSS. The cB98 form gives the best perfor¬ 
mance of all the functionals considered. It is noteworthy 
that in the mGGA class the mean errors reduce as more 
HF exchange is included in the functional - with cTPSS, 
cTPSS(h) and cB98 containing 0%, 10% and 19.85% HF 
exchange respectively. Suggesting that for this property 
that the treatment of exchange may be the dominant 
factor in the errors. Since the treatment of long-range 
exchange is not rectified in the transition from GGA to 
mGGA type functionals (and only partially corrected by a 
global admixture), then it may be that this factor far out 
weighs any improvements due to the inclusion of current 
effects. 

It is worth emphasizing that in the course of our inves¬ 
tigation we found that the implementation of the mGGA 
functionals including a generalized kinetic energy den¬ 
sity was straightforward. In particular we found that 
no special care was required with respect to numerics 
compared with standard functionals and that in practical 
use the functionals are robust and self-consistent calcula¬ 
tions using these functionals converge without significant 
difficulty. This sharply contrasts the behaviour for the 
VRG functionals as investigated in Refs. 16 and 17. 

The results for NMR shielding constants are presented 
in Figure 2. Here we see LDA is poor as expected. In 
addition KT3, which was designed for these properties, is 
the best performing functional - significantly improving 
over the standard PBE GGA functional and the B97-2 
hybrid functional. In this case we see the addition of the 
VRG(TP) correction to PBE and KT3 has little effect, 


100 


50 



FIG. 2. Box-whisker plot of errors in (C)DFT NMR shielding 
constants relative to the extrapolated CCSD(T) values of 
Ref. 43. All calculations use the aug-cc-pCVTZ basis set. See 
the text for further details. 

very slightly deteriorating the results. The mGGA results 
for cB98, cTPSS and cTPSS(h) are intermediate between 
PBE and KT3 - offering small systematic improvements 
over PBE. Again B98 produces the best results of the 
mGGA functionals. 

On the whole the quality of the mGGA results at mod¬ 
est field strengths may be regarded as disappointing. The 
overall errors suggest that mGGAs may offer modest im¬ 
provements over conventional GGA functionals such as 
PBE - though they cannot compete with GGAs tailored 
to specific properties. This is broadly consistent with 
findings by Bates and Furche for the calculation of exci¬ 
tation energies using cTPSS 13 . Since current corrections 
are known to be relatively small it is important that the 
underlying functional should be relatively accurate. For 
the mGGAs considered here there are known weaknesses 
(for example in the treatment of long-range exchange) 
that may obscure the effect of the current dependence. 
For the case of NMR shielding constants a more detailed 
analysis of the significance of current effects and how 
these interplay with errors in a range of density func¬ 
tionals is presented in Ref. 46. We will now examine 
how these functionals perform when the magnetic field 
becomes much higher and has a stronger effect on the 
electronic structure. 


B. The strong field regime: paramagnetic bonding 

One approach to explore whether or not the inclusion of 
current effects via the modified kinetic energy density of 
Eq. (13) is physically reasonable is to increase the strength 
of the magnetic field. Lange et. al. 44 have recently per¬ 
formed full configuration-interaction (FCI) calculations 
at high field that have uncovered a new mechanism for 
chemical bonding in the presence of a strong magnetic 
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TABLE III. Dissociation energies and equilibrium bond lengths 
for H 2 and rare gas dimers in a 1 a.u. magnetic field perpen¬ 
dicular to the internuclear axis. 



h 2 

Re 

He 2 

/ ao 
NeHe 

Ne 2 

h 2 

D e 

He 2 

/ mE h 
NeHe 

Ne 2 

HF 

2.708 

3.296 

3.543 

3.773 

2.340 

0.218 

0.482 

0.978 

LDA 

2.374 

2.550 

2.846 

3.062 

14.81 

8.231 

11.817 

19.730 

PBE 

2.514 

2.810 

3.080 

3.294 

5.985 

2.463 

3.938 

7.128 

KT3 

2.511 

2.852 

3.121 

3.342 

4.637 

2.661 

3.937 

6.527 

cB98 

2.640 

3.203 

3.472 

3.676 

1.328 

1.011 

1.514 

2.050 

cTPSS 

2.564 

2.864 

3.124 

3.346 

5.255 

1.307 

2.344 

4.577 

cTPSSli 

2.558 

2.879 

3.154 

3.379 

5.263 

1.245 

2.130 

3.990 

CCSD(T) 

2.578 

2.977 

3.248 

3.487 

4.554 

1.259 

2.217 

4.016 

FCF 

2.578 

2.975 

- 

- 

4.554 

1.271 

- 

- 


“ The He2 FCI calculations use the aug-cc-pVTZ basis set. 


field. This new bonding has been termed perpendicu¬ 
lar paramagnetic bonding and occurs at field strengths 
similar to those found on some white dwarf stars. Since 
this work Murdin et al. 4 ‘ have shown that phosphorus 
and selenium doped silicon semiconductors can produce 
a viable laboratory analogue of free hydrogen 4 ' and he¬ 
lium 48 in strong magnetic fields. The description of these 
types of systems via quantum chemistry will require less 
computationally demanding approaches - and CDFT is 
one strong candidate for the simulation of these systems. 

To investigate the performance of cB98, cTPSS and 
cTPSS(h) in strong magnetic fields potential energy pro¬ 
files were calculated for H 2 , He 2 , NeHe and Ne 2 - In 
particular, we consider the 3 £+(ler g l(7*) state of H 2 and 
the lowest 4 £+ states of He 2 , HeNe and Ne 2 - Each of 
these states is repulsive or weakly dispersion bound in 
the absence of a magnetic field but become more strongly 
bound when a field is applied. We note that only the 
3 £+(l<7glcr*) of H 2 is an overall ground state in the pres¬ 
ence of the field. These states were compared against 
results from accurate CCSD(T) potential energy curves 
calculated using a recent non-perturbative implementa¬ 
tion by Stopkowicz et al 29 in the London program. For 
comparison we have also generated similar profiles with 
standard LDA and GGA density functionals as well as 
with HF theory. The calculated potential energy curves 
for H 2 , He 2 , NeHe and Ne 2 are shown in Figures 3, 4, 5 
and 6, respectively. Equilibrium bond lengths and dis¬ 
sociation energies were determined numerically and are 
presented in Table III. 

The 3 £+(lcr ff l(7*) state of the H 2 molecule in a per¬ 
pendicular field was examined at the FCI level by Lange 
et al. A4: . The potential energy curves for this state are 
shown in Figure 3. HF strongly underbinds this state in 
comparison with the FCI data. In contrast LDA strongly 
overbinds, a tendency which is largely corrected by the 
PBE functional and further improved by the cTPSS and 
cTPSS(h) models. These trends are reflected in the equi¬ 
librium bond lengths and dissociation energies in Table III. 


Although not highly accurate the cTPSS and cTPSS(h) 
models give a reasonable qualitative description of the 
potential energy curve. The empirically parameterized 
KT3 functional is interesting because it gives simultane¬ 
ously a reasonable estimate of both the equilibrium bond 
length and dissociation energy. However, at intermediate 
separation an unphysical barrier is observed. For B98 
an even more pronounced barrier is present and the po¬ 
tential energy curve is generally even less accurate than 
Hartree-Fock theory. This may suggested that heavily 
parameterized functional forms, determined to perform 
well at zero field, may not be the best candidates for use 
in strong-field CDFT studies. 



FIG. 3. Potential energy curve for the 4 E+(l<7 9 lcq() state of 
the H 2 molecule in a magnetic Held of 1 a.u. perpendicular 
to the bonding axis for a variety of methods with the aug-cc- 
pCVTZ basis set. 

Examining the potential energy curves for He 2 in Fig¬ 
ure 4 we see that HF tends to under-bind with a bond 
length of 3.30 ao compared with the CCSD(T) value of 
2.98 a o- Similarly the HF dissociation energy of 0.218 
mEh is much smaller than the corresponding CCSD(T) 
value of 1.259 mEh- For this small system we were able 
to compare the CCSD(T) results with FCI values (cal¬ 
culated in the slightly smaller aug-cc-pVTZ basis), as 
expected the agreement is excellent - the corresponding 
potential energy curves are essentially indistinguishable 
on the scale of Figure 4. For LDA we see a strong ten¬ 
dency to overbind giving much too short R e values and 
much too large estimates of D e . The GGA functionals 
PBE and KT3 show considerable improvement over LDA, 
however, they still strongly overbind. The improvement 
for the mGGA functionals is striking - in particular TPSS 
and TPSS(h) give a good qualitative description of the 
potential energy curve. The corresponding R e and D e 
values indicate that there still remains a tendency towards 
over binding but this is greatly reduced. 

The cB98 functional tends to show more significant 
under binding. Here we note that the arguments used in 
the construction of cB98 and cTPSS are rather different. 
In particular, cB98 is an empirically parameterized func¬ 
tional (see the appendix), whereas cTPSS is constructed 
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FIG. 4. Potential energy curve for the lowest 1 E^' state of He 2 
in a magnetic field of 1 a.u. perpendicular to the bonding axis 
for a variety of methods with the aug-cc-pCVTZ basis set. 

based on the satisfaction of known exact conditions. In 
this work we have used the parameters determined in 
Ref. 27 to define the cB98 form. These parameters were 
determined at zero field and from post-LDA calculations - 
as a result they may not be optimal for fully self-consistent 
calculations in the presence of a magnetic field. On the 
other hand the cTPSS functional is designed to satisfy se¬ 
lected constraints at zero field and it could be argued that 
in the presence of a field both the B98 and TPSS based 
functionals are open to further optimization, though this 
is beyond the scope of the present work. 

The stability of the mGGA functionals is particularly 
evident in the the strong field regime when one compares 
the present results for He 2 with those for the VRG-based 
estimates in Figure 7 of Ref. 17. The VRG approaches led 
to very difficult SCF convergence and complex potential 
energy curves with a strong unphysical over binding. The 
mGGAs considered here are un-problematic in practical 
application and yield results surprisingly close to the 
CCSD(T) estimates. 

Similar qualitative trends are observed for the NeHe and 
NeNe dimers in Figures 5 and 6. Again LDA and GGA 
functionals are not sufficiently accurate for practical use 
and the mGGA functionals provide a large improvement. 
The cTPSS and cTPSS(h) results remain impressive - 
with cTPSS (h) being consistently slightly more accurate 
than cTPSS. This trend is reflected in both the potential 
energy curves and the R e and D e values in Table III. 


V. INTERPRETATION OF PARAMAGNETIC BONDING 
IN THE KS-CDFT FRAMEWORK 

The mGGA CDFT functionals offer a computationally 
cheap correlated method for the examination of the ex¬ 
otic bonding mechanisms observed in a strong magnetic 
field. In many areas of chemistry the nature of bonding, 
chemical reactions, spectra, and properties of molecular 
species are interpreted qualitatively in terms of orbital 



FIG. 5. Potential energy curve for the lowest 1 E;J" state of the 
NeHe dimer in a magnetic field of 1 a.u. perpendicular to the 
bonding axis for a variety of methods with the aug-cc-pCVTZ 
basis set. 



FIG. 6. Potential energy curve for the lowest 1 Ej' state of Ne 2 
in a magnetic field of 1 a.u. perpendicular to the bonding axis 
for a variety of methods with the aug-cc-pCVTZ basis set. 

interactions. We now consider the extent to which in¬ 
formation from KS-CDFT calculations can aid in simple 
interpretation of the perpendicular paramagnetic bonding 
interactions. 

We begin by considering a molecular orbital analysis 
of the perpendicular paramagnetic bonding. KS-CDFT 
calculations provide a simple set of canonical molecular 
orbitals, which can be used to construct the electronic 
density via 

p( r ) = l^( r M( r )l 2 = J2 U} 'y( r ') D 'v c w c( r ) ( 16 ) 

i 7C 

Here we note that the occupied KS orbitals <^(r) can 
be complex in the presence of a magnetic field. In the 
second equality the density is expressed in terms of the 
one-particle density matrix D 7 £ and the basis functions 
w 7 (r). We will commence by considering how the molec¬ 
ular orbital energies associated with H 2 and the rare gas 
dimers change upon application of a magnetic field as 




































the perpendicular paramagnetic bonding in Section V A 
evolves. Since the orbitals themselves can be complex in 
the presence of a field we then proceed in Section V B to 
analyze the bonding in terms of the changes in electronic 
density of Eq. (16) as a function of field, which is naturally 
a real observable quantity. 



w 

£ 


H Is 

/ 1.301 

\ His 

A. KS-CDFT molecular orbital analysis 

- 10.0 

0.0 

°g 

-1.325 

0.0 



KS molecular orbitals have been widely used as an 
interpretive aid in chemical applications throughout the 
literature. The KS orbitals are defined to minimize the 
non-interacting kinetic energy and yield the physical elec¬ 
tronic density via Eq. (16). They also have appealing 
properties; for example the highest occupied MO energy 
is minus the first ionization potential (IP) 49,50 and the re¬ 
maining orbital energies can be interpreted as Koopman’s 
type approximations to higher IPs 51 . The extent to which 
these properties hold for general practical approximations 
has been a subject of debate in the literature 52,53 , as 
has the interpretation of KS virtual orbitals 54 owing to 
the role of the integer discontinuity 49 , which is missing 
from common approximations. However, from a practical 
standpoint it is widely accepted that interpretations based 
on occupied KS orbitals (to which we limit the following 
discussion) are theoretically justified and their utility has 
been borne out in many practical applications. 

We now consider how the KS orbital energies change 
upon application of a perpendicular magnetic field of 1 a.u. 
Given that the cTPSS based models seem to be the most 
reliable of those studied in the present work we consider 
how the orbital energies from this functional change in 
Figures 7 and 8 for H 2 and He 2 , respectively. For the H 2 
molecule in the 3 E+(ler g l<T*) state we consider the energy 
of the occupied a g and a* orbitals change upon bonding. 
In particular, we plot orbital energies in the absence of 
a field and in a perpendicular field of 1 a.u. relative to 
those of the atomic orbitals in the same field. We have 
plotted the orbital energies at an internuclear separation 
R e = 2.564 a.u. consistent with the cTPSS equilibrium 
bond length in the presence of a field. 

We see that in the absence of a field the singly occupied 
<jg orbital is stabilized by 1.33 mEj,, whilst the singly 
occupied a* is destabilized by 1.30 mEh relative to the 
Is hydrogen orbitals. This is consistent with a net bond 
order of zero and a repulsive profile for the corresponding 
potential energy curve. In the perpendicular field of 1 a.u. 
we see that, relative to the hydrogen Is orbital in the same 
field, the a g orbital is stabilized by 63.21 mEh, whilst the 
a* orbital is destabilized by 34.97 mEh- This greater 
stabilization of the a g orbital leads to a net bonding 
interaction, consistent with the analysis of Ref. 44. This 
illustrates that the KS orbitals can be a useful tool in 
rationalizing this exotic bonding phenomenon. 

A similar analysis can be carried out for the lowest 1 E+ 
state of He 2 and is presented in Figure 8. In this plot 
we have separated the spin down and spin up orbitals 


<V 


B, 1 


34.971 


B , 0 


-63.207 

FIG. 7. Relative orbital energies for the occupied a g and 
a* CDFT molecular orbitals for H 2 relative to the atomic 
Is orbitals, with (blue) and without (red) a field of 1.0 a.u. 
perpendicular to the interatomic axis. 


and defined their energies relative to atomic orbitals of 
the same spin in the same field. Defined in this way 
the offset between the orbitals of different spin due to 
the Zeeman interaction is removed from the plot. Again 
the energies correspond to the cTPSS He 2 equilibrium 
internuclear separation R e = 2.864 a.u. in the presence of 
a perpendicular field. In the absence of a field we see that 
again the relative stabilization of er g and destabilization of 
ct* are approximately compensatory, whilst in the presence 
of a field the <r g orbital is more stabilized than the a* 
orbital is destabilized. It is also clear that the extra 
stabilization of the cr g orbital of ~ 8 mEj, is considerably 
less than the ~ 28 mEh observed for H 2 . This is consistent 
with the strength of binding exhibited for these species 
in Figures 3 and 4 as well as in Table III. 

Similar orbital energy diagrams may be constructed 
for the HeNe and NeNe systems, however, they become 
significantly more complex due to large differences be¬ 
tween the orbital energies in HeNe and the splitting of the 
p-orbitals in NeNe. We therefore consider an alternative 
visualization of the bonding effects in these systems based 
on the charge and (physical) current density differences. 


B. Electron density analysis 

For each of the species He 2 , HeNe and Ne 2 , we have per¬ 
formed calculations of the electronic density and current 
density in the presence of varying perpendicular magnetic 
fields (B±) using the cTPSS functional at the correspond- 
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He Is 1 
0.0 


°u" B ± 0 

2.428 


B , 0 


2.428 


-2.473 


-2.473 


FIG. 8. Relative orbital energies for the occupied o g and 
<r* CDFT molecular orbitals for He 2 relative to the atomic 
Is orbitals, with (blue) and without (red) a held of 1.0 a.u. 
perpendicular to the interatomic axis. 


ing equilibrium geometries. Here we consider the density 
change Aps x (r) for each system relative to the isolated 
atoms in the same field, 


A PB± (r) = p^(r) - £ Pb ± ( r ) (17) 

2=1 

This difference density allows for a visualization of the 
nature of the paramagnetic bonding in these systems. In 
a similar manner one can consider the (gauge invariant) 
physical current density difference Aj 

N atoms 

Aj S ,(r)=j° imer (r)- £ j‘ Bx (r) (18) 
2=1 

In Figure 9 we present plots of the density differences in 
Eqs. (17) and (18) for each of the species with B± = 1.0 
a.u. The shading of the contours represents the buildup 
(red) or depletion of the charge density (blue), relative 
to two non-interacting atoms in the same field. In all 
three cases there is a clear build up of density between 
the atoms consistent with bonding. The charge density 
difference is elongated along the held, above and below 
the plane of the plots. The streamlines show the vector 
Held associated with the current density difference. Para- 
tropic circulations are clearly visible over the centre of 
the bonds where charge density accumulates, and diat- 
ropic circulations are visible in regions where the charge 
density is depleted. We have confirmed that for higher 
Helds the alignment between para- / dia-tropic circula¬ 


tions and charge accumulation / depletion becomes more 
pronounced. 

This picture of the perpendicular paramagnetic bond¬ 
ing suggests that as the charge density is elongated along 
the field the constituent atoms may approach one another 
more closely. As they do so they experience a greater 
nuclear-nuclear repulsion, which may be screened by a 
rearrangement of the charge density towards the bond 
centre. This rearrangement is accompanied by consistent 
para- and dia-tropic current circulations. The cTPSS 
functional used in this work provides a simple, computa¬ 
tionally cheap, route to perform analysis of the bonding 
encountered in the strong field regime. 


VI. CONCLUSION 

In this work we have implemented a previously 
detailed 23,24,34 modification to the kinetic energy density 
term in mGGA functionals to perform non-perturbative 
cDFT calculations in a stable, gauge invariant manner 
using the London program 18 . The modified mGGA 
functionals cB98, cTPSS and cTPSS(h) show a level of 
accuracy in predicting weak field magnetic properties that 
is competitive with existing GGA functionals without any 
additional fitting. The functionals cTPSS and cTPSS(h) 
show excellent prediction of perpendicular paramagnetic 
bonding behaviour, suggesting that the modification of 
the kinetic energy density is a viable route for the incor¬ 
poration of current effects in standard mGGA density- 
functionals. In contrast to vorticity dependent forms 
previously studied 16,1 ' the functionals exhibited excellent 
numerical stability in the finite field setting without the 
need for delicate numerical regularizations. 

Whilst the mGGA results show considerable promise 
in the high field regime, their performance in the weak 
field regime is perhaps disappointing - leading to only 
modest improvements of conventional GGA forms. To 
some extent this may be due to the fact that mGGAs 
contain many of the shortcomings associated with GGA 
forms. For example, the functionals still have potentials 
with incorrect asymptotic behaviour - particularly for the 
exchange contribution. Since the current effects at weak 
field are not dominant but rather add small corrections to 
the predominantly Coulombic exchange and correlation 
interactions then it may be necessary to further improve 
the underlying functional forms before the true impact of 
the current terms can be assessed in this regime. 

We expect that this approach to include current depen¬ 
dence should play a central role in the future development 
of new CDFT functionals and many avenues are open for 
development. Obvious possibilities include the generaliza¬ 
tion of range-separated mGGA functionals 55,56 to obtain 
a better balance of errors between exchange, correlation 
and current contributions and the re-parameterization of 
functionals either empirically or via the consideration of 
alternative exact conditions in their construction. In the 
latter category we note that the presence of a magnetic 
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FIG. 9. ApB ± (r) (coloured contours) and Ajs^(r) (streamlines) for the rare gas dimers He 2 (left), HeNe (middle) and Ne 2 
(right) in a B± = 1.0 a.u. magnetic Held. The internuclear axis is aligned with the z-axis and the plots are show in the yz -plane 
intersecting the atomic positions. 


field causes compression of the electronic density in two 
dimensions perpendicular to the field and elongation along 
it. As a result non-uniform coordinate scaling relations 
may be one example of conditions that may provide a 
powerful tool in further development. The London pro¬ 
gram 18 provides a powerful platform for this work since 
CDFT approaches can be calibrated against accurate ab 
initio data for a range of field strengths, where experi¬ 
mental data may be scarce. In the future we hope that 
these relatively inexpensive CDFT approaches may then 
be applied to the study of larger systems such as those in 
Refs. 47 and 48. 


Appendix A: meta-GGAs in CDFT 

We present the key working equations defining the B98 
and TPSS functionals, indicating how the the modified 
T a of Eq. (13) enters each functional. 


where the exchange hole curvature for the uniform electron 
gas takes the simple form 

Q^G = _l (67r 2 )2 /3 p 5/3 (A2) 

and Q c , is defined in Eq. (10). This inhomogeneity fac¬ 
tor, q a , controls the enhancement or attenuation of the 
exchange and correlation energy over the uniform gas 
values. 

The exchange component of the functional takes the 
form 

E x ,a = J e™ G (p a )g x , a (q a )dr (A3) 

and the opposite- and same-spin components of the cor¬ 
relation energy take the forms 

E c ,otp = / e c,a/3 {Pon P/3)Sc,a/3(<7avg)dr (A4) 


1. cB98 


E, 


= / e^{p a )fa gc,aa(qa)dr. 


(AS) 


The B98 functional 2 ' has a general construction that 
is similar to the popular B97 functional form 57 , however, 
instead of using reduced spin-density gradients it makes 
use of an inhomogeneity parameter q a 


The functions g x , a (qa), 3c,a/3(<?av g ) and g c ^ a (q a ) are di¬ 
mensionless inhomogeneity correction factors depending 
on q a and <? avg = \(q a + qp). In addition the same-spin 
correlation energy contains a self-correlation correction 
(SCC) factor 


( q * - Qr G ) 

IQ GEG I 


f! cc 


1 (Vp ff ) 2 


qcr = 


(Al) 


4 Pc r 


(A6) 
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This factor varies between 0 and 1 and vanishes in one- 
orbital regions, ensuring the functional is self-correlation 
free. For convenience in deriving fitted forms for the 
inhomogeneity correction factors g , Becke proposed the 
following transformation to a finite interval, 


79 

y/l + 7 V 


(A7) 


where 7 is a parameter to be determined and separate 
transformations are carried out for the exchange, opposite- 
spin correlation and like-spin correlation respectively us¬ 
ing the appropriate definitions of q as in Eqs (A3)-(A4). 
Based on calculations of atomic exchange-correlation en¬ 
ergies Becke proposed the values 7 Xi(T = 0.11, 7 c ,«/3 = 0.14 
and 7 c,o-o- = 0.16. The inhomogeneity corrections g were 
then determined by fitting a power series expansion of 
the form 


<? = £ CiW 1 (A 8 ) 

2=0 


with to = 2 chosen to prevent unphysical over-fitting. 
This fitting was carried out using the G2 thermochemical 
dataset and basis set free, post-LSDA, calculations. The 
optimal parameters can be found in Ref. 27. 

Here we use this parameterization directly but note that 
in future these parameters could be re-optimized based 
on self-consistent data. In addition an amount of Hartree- 
Fock exchange is included with weight c x = 0.1985. The 
resulting B98 functional therefore possesses a high degree 
of non-locality and may be classified as a hybrid mGGA 
functional with dependence not only on r but also on the 
laplacian of the density V 2 p. 

The original definition of B98 utilized the zero-field 
exchange hole curvature of Eq. (10) in its definition. How¬ 
ever, it was noted 2 ' that this form can be readily extended 
to include current effects via Eq. (11) and it is this avenue 
that we explore in the present work. Unless otherwise 
stated we employ this modified form throughout and 
denote it as cB98. 


2. cTPSS 


One of the most widely used meta-GGA functionals is 
due to Tao, Perdew, Staroverov and Scuseria (TPSS) 26 . 
This functional is designed to satisfy exact constraints 
without empirical parameters and as such is an interesting 
candidate to study in the context of generalization to 
finite magnetic field strengths where much less is known 
about the performance of approximate functionals. In 
this functional the ratio 

z = 2t w /t , r = ^r CT , r w = l!^L (A9) 

8 0 

rr ' 


plays a key role as a dimensionless inhomogeneity param¬ 
eter, along with 


IVPI 2 2 

4(3tt 2 ) 2 /V/3 


(A10) 


Note that throughout this work we use the definition of 
r in Eq. (9), which does not include the factor of 1/2 
commonly employed. The exchange functional then takes 
the form 


E x [p\ = J p£ BEG (p)F x {p,z) (All) 

where the precise details of the form chosen for the en¬ 
hancement factor F x (p, z ) can be found in Ref. 26. The 
correlation energy takes the form 


E c [PonPp, Vp a , S7pp,r] = 

J p£ T c evPKZB (p a , p 0 ,S7 p ai S7 p p ,t) 

x [l+ de™ vPKZB (/ 9 a,P/ 3 ,Vp a ,Vp^,r)(r w /r) 3 ] dr 

(A12) 


where d = 2.8 hartree -1 . 

Using the replacement in Eq. (13) leads to modifications 
in the exchange contribution via z in Eq. (A9) and in the 
correlation energy as shown in Eq. (A12). This modified 
form is noted cTPSS and is consistent with that used in 
the response implementation of Ref. 13. 
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